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We  solve  the  inviscid  Euler  equations  in  complicated  geometries  using  a  Cartesian 
grid.  This  requires  solid  wall  botmdeiry  conditions  in  the  irregular  grid  cells  near  the 
boundary.  Since  these  cells  may  be  orders  of  magnitude  smaller  than  the  regular  grid 
cells,  stability  is  a  primary  concern.  We  present  a  new  approach  to  this  problem  and 
illustrate  its  use.  ( 
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1  Introduction 


In  previous  work  [1],  [6],  [7],  we  have  described  a  Cartesian  grid  method  for  the 
inviscid  Euler  equations  in  arbitrary  geometries.  There  are  many  advantages  to  be 
gained  from  this  approach.  Grid  generation  is  simplified,  since  we  avoid  the  use  of 
(possibly  multiblock)  body-fitted  grids,  and  we  can  use  high  resolution,  highly  efficient 
solvers  on  regular  grids  over  the  bulk  of  the  domain.  This  has  led  to  renewed  interest 
in  Cartesian  grids  in  recent  years,  e.g.,  [3],  [10].  One  of  the  difficulties  with  Cartesian 
grids  is  that  they  give  insufficient  resolution  in  certain  regions  such  as  leading  edges. 
This  can  now  be  overcome  by  Cartesian  adaptive  mesh  refinement  [1],  [2]. 

The  principal  remaining  difficulty  in  this  approach  is  due  to  the  essentially  ar¬ 
bitrary  way  that  a  Cartesian  grid  intersects  the  boundaries  of  the  computational 
domain.  In  particular,  a  solid  wall  boundary  cutting  through  the  grid  creates  irreg¬ 
ular  cells  that  may  be  orders  of  magnitude  smaller  than  the  regular  cells  away  from 
the  boundary.  For  these  irregular  cells,  special  difference  equations  are  needed  that 
maintain  stability  and  accuracy,  and  satisfy  the  solid  wall  boundary  conditions  of  no 
normal  flow. 

In  this  work,  we  present  an  improved  method  for  the  small  boundary  cells.  We 
use  an  explicit,  finite  volume  formulation  that  computes  fluxes  at  cell  edges  on  the 
regular  part  of  the  domain.  We  would  like  to  define  fluxes  at  the  edges  of  the  irregular 
cells  in  such  a  way  that  the  method  is  stable  even  with  very  small  boundary  cells, 
using  a  time  step  based  on  the  regular  grid  cells  away  from  the  boundary.  The  CFL 
condition  requires  that  the  numerical  method  allow  information  to  propagate  at  least 
as  quickly  as  the  underlying  differential  equation.  In  the  present  context  this  means 
that  we  must  define  fluxes  at  the  sides  of  our  irregular  cells  based  on  more  than  just 
the  neighboring  cell  values. 

In  our  previous  work,  we  have  used  a  wave  propagation  approach  in  defining  these 
fluxes.  Here  we  propose  an  alternative  method  that  has  some  advantages  over  the 
wave  propagation  approach.  In  particular,  the  wave  propagation  method  is  subject 
to  intermittent  instabilities  due  to  two-dimensional  effects  that  are  not  clearly  under¬ 
stood.  The  new  method  has  a  cancellation  property  in  two  dimensions  that  appears  to 
give  better  stability  properties.  Moreover,  the  computational  geometry  is  simphfied 
in  the  new  approach.  The  fluxes  are  defined  in  terms  of  weighted  averages  of  nearby 
cell  values.  These  weights  may  be  calculated  as  a  preprocessing  step  on  any  fixed  grid 
and  need  not  be  repeatedly  calculated.  In  the  previous  approach  the  weights  depend 
on  the  flow  variables  and  a  certain  amount  of  computational  geometry  was  required 
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near  the  boundary  in  every  time  step. 

We  consider  the  inviscid  Euler  equations  of  gas  dynamics  in  two  space  dimensions, 


ut  +  /(^i)x  +  ^(u)v  =  0 


(1.1) 


where 


p 

pui 

pU2 

u  = 

pUi 

pU2 

.  /(l^)  = 

pu\+p 

PUIU2 

,  giu)  = 

PU1U2 
pu\  -f  p 

pE 

.  ui{pE+p)  . 

.  U2{pE-]-p)  . 

Here  {ui,U2)  represents  the  velocity,  E  is  the  total  energy  per  unit  mass,  and  p  is  the 
pressure,  which  is  related  to  the  other  variables  by  the  equation  of  state.  We  assume 
a  7-law  gas,  so  that 

P  =  (7  -  +  Wa)).  (1.3) 

At  a  solid  wall  boundary  we  require  that  the  component  of  velocity  normal  to  the 
wall  be  zero. 

In  one  space  dimension  the  system  reduces  to 

^^^+fi^t)x  =  0  (1.4) 

where  u  =  [p,pv,pE)  and  f{u)  =  (pv,pv^  +p,v{pE  +  p)),  with  v  =  Ui  the  velocity. 
The  boundary  conditions  become  v  =  0  at  a  solid  wall. 


2  A  one- dimensional  example 

In  order  to  illustrate  this  approach  we  begin  with  a  one-dimensional  model  problem, 
the  one- dimensional  Euler  equations  for  x  >  0  with  a  solid  wall  at  s  =  0.  We  take  a 
grid  with  cell  interfaces  at  the  points 


xo  =  0 
Xi  =  h' 

Xj  =  h'  +  jh  for  j  =  2,  3 . 


Here  h  is  a.  uniform  grid  spacing  and  h'  <  h.  the  grid  is  uniform  except  for  one  small 
cell  near  the  boundary  (see  Figure  1).  We  use  a  conservative  method  in  the  form 


- 


k 


=  3=0,1, 


(2.5) 
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Figure  1:  One  dimensional  grid  with  one  irregular  cell  adjacent  to  the  wall. 

Here  hj  is  the  width  of  the  jth  cell,  so  in  our  case  we  have  ho  =  h'  and  hj  =  h  for 

j  >  0. 

For  simplicity  we  restrict  our  attention  to  Godunov’s  method  on  the  regular  por¬ 
tion  of  the  grid,  although  the  ideas  we  propose  can  be  extended  to  higher  order 
methods  cis  well.  In  Godimov’s  method  we  take 

(2-6) 

where  represents  the  solution  to  the  Riemann  problem  with  left  and  right 

states  and  u^,  evaluated  along  i/t  =  0.  Although  a  rigorous  stability  proof  is  not 
available  for  systems  of  equations,  in  practice  this  method  is  always  stable  provided 
the  CFL  condition 


is  satisfied,  where  is  the  maximum  wave  speed.  We  will  assume  that  our  time 
step  k  is  chosen  so  that  the  condition  (2.7)  is  satisfied  relative  to  the  uniform  h.  We 
will  use  the  flux  (2.6)  for  j  =  2,  3,  . ..,  i.e.,  at  all  interfaces  where  the  cell  on  both 
sides  is  regular.  Our  task  is  to  define  fluxes  for  j  =  0,  1  so  that  we  maintain 
stability  (emd  accuracy)  with  this  time  step  even  if  h'  <C  h. 

First  suppose  h'  =  h.  Then  we  can  use  the  Godunov  flux  (2.6)  zdso  at  j  =  1. 
At  the  wall  we  use  the  well-known  observation  that  the  solution  to  the  boundary 
value  problem  can  be  obtained  by  ignoring  the  wall  wd  extending  the  computational 
domain  to  the  whole  line  — oo  <  i  <  oo,  if  we  take  data  uo{x)  for  x  <  0  equal  to 

p(x,0)  =  p(-x,0) 

v(x,0)  =  — v(— x,0)  for  X  <  0 

p(x,0)  =  p(-x,0). 

We  denote  this  “reflection”  of  the  data  (in  which  the  velocity  is  negated)  by  the 
operator  71,  so  that  for  shorthand  we  can  write 

u(x,  0)  =  7^(u(-x,0))  for  X  <  0. 
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With  this  extended  data,  the  solution  continues  to  satisfy  u(x,t)  =  7Z(u(—x,t))  also 
for  t  >  0  and  in  particular  the  boundary  condition  u(0,t)  =  0  is  automatically 
satisfied.  This  suggests  that  we  obtain  a  flux  at  the  wall  by  solving  a  Riemann 
problem  with  left  and  right  states 

=  7Z(Uo),  =  Uo 

in  each  time  step  to  obtain 


Fo  =  fiu*{n{Uo),Uo). 


(For  brevity  we  will  leave  off  the  superscript  n  in  general.)  Note  that  the  density 
and  energy  components  of  this  flux  will  be  zero  since  the  velocity  component  of  u*  is 
zero  at  the  wall.  There  will  only  be  a  momentum  flux  at  the  wall  due  to  the  pressure 
there,  as  expected  physically. 

li  h'  <  h  we  could  attempt  to  use  this  same  formula  to  define  Fq  but  we  would 
find  that  it  is  unstable  unless  the  CFL  condition 


kXr 


h> 


<  1 


(2.8) 


is  satisfied.  This  will  place  an  unreasonable  restriction  on  h. 

This  instability  is  caused  by  the  fact  that  the  boundary  flux  Fq  is  based  on  the  data 
Uq  alone.  If  the  CFL  condition  (2.8)  is  satisfied,  then  it  is  only  this  data  that  affects 
the  flux  at  the  wall  over  the  time  step.  However,  when  (2.8)  is  violated  the  value 
should  also  affect  the  flux  at  the  wall,  and  ignoring  this  effect  leads  to  instability. 

In  a  “large  time  step”  approach  we  increase  the  stencil  of  the  method,  meaning  we 
allow  more  data  points  to  come  into  the  computation  of  each  flux,  and  hence  retain 
stability.  One  way  to  achieve  this  is  by  a  wave  propagation  approach.  The  solution  of 
the  Riemann  problem  at  each  cell  interface  consists  of  three  waves  propagating  away 
from  the  interface.  If  (2.8)  is  satisfied  then  these  waves  remain  in  the  cells  bordering 
the  interface  during  the  entire  time  step  and  hence  affect  the  solution  only  in  these 
cells.  If  (2.8)  is  violated  then  the  waves  may  affect  cells  further  away.  Implementing 
Godunov’s  method  in  terms  of  this  wave  prppagation  approach,  allowing  waves  to 
affect  more  than  just  the  adjacent  cell,  gives  a  large  time  step  generalization  that 
remains  stable  for  much  larger  time  8teps[5].  In  the  present  context  this  allows  us 
to  reduce  h'  without  reducing  the  time  step  k.  Waves  from  the  boundary  Riemann 
problem  cross  the  interface  at  xj  and  affect  Ui  as  well  as  Uq.  Waves  from  the  interface 
at  Xi  may  reach  the  boundary.  These  waves  reflect  off  the  boundary  and  the  reflected 
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«vave  affects  the  value  Uq  and  perhaps  also  U\  if  the  reflected  wave  reaches  the  cell 
interface  at  xi  during  the  time  step. 

A  more  detailed  description  of  this  procedure  may  be  found  in  [7].  A  natural 
extension  to  two  space  dimensions  gives  one  method  to  deal  with  small  cells  near 
the  boundary,  as  described  in  [1],  [6],  [7].  In  one  dimension  this  works  very  well 
but  in  two  dimensions  occasional  stability  problems  have  still  been  observed  due  to 
multidimensional  effects. 

The  nev'  approach.  Our  new  approach  to  the  small  cell  problem  can  also  be 
illustrated  with  the  one- dimensional  problem  described  above.  We  again  use  the 

method  (2.5)  with  Godunov  fluxes  (2.6)  for  j  =  2,  3,  _  For  j  =  0  and  ji  =  1  we 

define  fluxes  in  a  similar  manner  but  with  a  different  choice  of  states  and  in  the 
Riemann  solver.  Recall  that  in  a  naive  attempt  to  use  Godunov’s  method  regardless 
of  the  size  h'  of  the  small  cell  we  would  take  left  and  right  states 

=  Tl(Uo)  <  =  Uo  (2.9) 

uf  =  I/o  u^  =  Ui.  (2.10) 

To  maintain  stability  when  h'  is  small,  we  need  to  allow  data  from  additional  grid 
cells  to  affect  the  left  and  right  states  at  each  of  these  interfaces.  Recall  that  the 
method  is  assumed  to  be  stable  with  our  choice  of  k  and  h  on  the  regular  portion  of 
the  grid.  This  suggests  that  we  should  define  by  taking  the  average  value  of  U 
over  an  interval  of  length  h  to  the  left  of  the  interface  Xj  and  define  u"  by  taking  the 
average  value  of  U  over  an  interval  of  length  h  to  the  right  of  Xj. 

For  example,  at  xq  =  0  (the  wall)  we  set 

u^  =  hh'Uo  +  {h-h')U{)  (2.11) 

ti 

If  we  view  the  grid  values  as  defining  a  piecewise  constant  function  with  values  Uj 
in  the  jth  cell,  then  (2.11)  is  the  average  value  of  this  function  over  the  interval 
0  <  X  <  h.  Note  that  ii  h'  =  h  (the  grid  is  Completely  regular)  then  (2.11)  reduces 
to  Uq  =  Uq  as  expected  for  Godunov’s  method.  Recall  that  in  Godunov’s  method 
we  take  Uq  =  7l(Uo)  =  7^('Uo)  impose  the  boundary  condition  u(0,t)  =  0.  This 
suggests  that  more  generally  we  take 

<  =  (2.12) 

where  itj  is  defined  by  (2.11).  We  then  use  the  Godunov  flux 

i^o  =  /(u*(«))  (2.13) 
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as  the  flux  at  the  wall.  Using  (2.12)  guaxantees  that  there  will  be  no  flux  of  mass  or 
energy  through  the  wall  and  hence  that  the  method  is  conservative. 

To  define  the  left  and  right  states  at  xi  we  again  construct  intervals  of  length  h 
to  either  side  of  this  point  and  average  the  piecewise  constant  function  defined  by  U 
over  these  intervals.  To  the  right  of  xi  lies  a  regular  cell  of  length  h  and  so 

uf  =  Ui.  (2.14) 

To  the  left  of  Xi  an  interval  of  length  h  extends  beyond  the  wall  (assuming  h'  <  h). 
Beyond  the  wall  we  assume  that  U  takes  the  value  Uq  given  by  (2.12).  A  weighted 
average  of  this  value  and  Uq  gives  uf: 

=  hh'Uo  +  ih-  h')u^o)-  (2-15) 

h 

The  flux  fi  is  then  defined  by 

F,  =  (2.16) 

Again,  \i  h'  =  h  this  reduces  to  the  standard  Godunov  flux. 

This  method  remains  stable  even  when  h'  h.  To  see  why  this  should  be  so, 
consider  the  formula  (2.5)  for  j  =  0  where  hj  =  k'.  It  is  the  division  by  h'  that 
may  cause  stability  problems  unless  the  fluxes  Fq  and  Fi  themselves  agree  to  0{h') 
as  h'  -*  0.  The  Godunov  fluxes  based  on  (2.9),  (2.10)  do  not  have  this  property. 
However,  our  proposed  fluxes  (2.13)  and  (2.16)  do  have  this  property,  since  inspection 
of  the  formulas  (2.11),  (2.12),  (2.14),  and  (2.15)  shows  that  ttf  =  uj  +  0{h')  jind 
uj  =  uf  +  0{h')  as  h'  — »  0.  Since  the  fiiix  function  f{u*{u^,u^))  is  a  Lipschitz 
continuous  function  of  and  u”,  it  follows  that  Fi  —  Fq  =  0{h')  as  /i'  — »  0  and  there 
is  at  least  a  chance  that  the  method  remains  stable  for  arbitrary  h'  <C  h.  Numerical 
experiments  show  that  this  is  indeed  the  case  (although  it  is  possible  to  contrive 
examples,  such  as  a  strong  rarefaction  wave  originating  at  this  irregularity,  where  the 
results  are  not  very  accurate). 

3  Boundary  conditions  in  two  dimensions 

Turning  now  to  the  two-dimensional  problem,  we  will  give  a  brief  description  of  how 
the  idea  described  above  extends  to  handle  the  small  cell  problem. 

Consider  the  portion  of  the  boundary  shown  in  Figure  2a  and  a  typical  boundary 
cell  (i,  j).  The  formula  for  updating  the  value  Uij  is  the  two-dimensionzd  analog  of 
(2.5), 
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Figure  2:  (a)  The  Cartesicin  grid  near  the  boundeiry.  (b)  Blow  up  of  cell  (ij)  showing 
the  location  of  fluxes. 

The  fluxes  F,  G,  and  H  represent  flux  per  unit  time  through  the  corresponding  side 
of  the  grid  cell  (see  Figure  2b)  and  Aij  is  the  area  of  the  cell.  If  any  of  the  sides  are 
missing,  the  corresponding  flux  is  zero. 

On  regul<ir  grid  cells,  Hij  =  0  and  the  fluxes  F  and  G  might  be  defined  by  an 
extension  of  the  Godunov  method,  setting 

Fj  =  Uii)),  Gii  =  hg{y;{Ui^.u  Ua)).  (3.18) 

Here  u*  represents  the  solution  to  the  appropriate  one- dimensional  Riemann  problem 
in  the  i  or  y  direction.  Note  that  the  fluxes  include  the  factor  h,  the  length  of  each 
side,  to  give  a  flux  per  unit  time  across  the  side. 

It  is  the  denominator  Aij  in  (3.17)  that  causes  trouble  when  the  cell  is  very  small. 
We  again  assume  the  method  is  stable  on  the  regulw  portion  of  the  grid,  where 
Aij  =  A*.  To  maintain  stability  we  need  to  insure  that  our  formulas  for  the  fluxes 
cause  the  total  flux  (the  sum  in  brackets  in  (3.17))  to  cancel  to  0{Aij)  as  Ajj  — >  0. 
This  is  only  possible  if  the  fluxes  are  computed  via  formulas  that  involve  more  than 
just  the  two  cells  bordering  the  cell  side.  We  take  an  approach  emalogous  to  what  we 
described  above  in  one  dimension. 

Boundary  fluxes.  We  begin  by  considering  the  boundary  segment,  where  we 
must  compute  the  flux  Hij.  In  two  dimensions  the  solid  wall  boundary  condition 
requires  that  the  normal  velocity  at  the  wall  be  equal  to  zero.  If  we  have  some  value 
■uj“  representing  the  value  of  u  just  inside  the  wall,  then  we  can  obtain  the  flux  Hij 
by  solving  a  one- dimensional  Riemann  problem  in  the  direction  normeil  to  the  wall, 
with  left  and  right  states 


*3 


u-  =  u‘“. 

*3 
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Figure  3:  The  inbox  and  outbox  constructed  from  the  boundary  segment  of  cell  (t, 
and  the  inbox  for  two  neighboring  cells. 

The  reflection  operator  7?.  is  now  defined  by  negating  the  normal  velocity  component 
while  leaving  the  tangential  velocity  component  along  with  the  density  and  pressure 
unchanged.  The  resulting  Godunov  flux  is  used  for  Hij. 

We  obtain  by  a  procedure  analogous  to  the  one-dimensionaJ  example.  We 
construct  a  box  extending  a  distance  h  away  from  the  wall  as  shown  in  Figure  3.  The 
box  extending  into  the  computational  domain  is  Ceilled  inbox(t,_;).  The  mirror  image 
box  outside  the  domain  is  called  outbox(i,j).  We  obtain  the  value  u‘>  by  viewing 
the  given  data  U  as  defining  a  piecewise  constant  function,  constant  in  each  grid  ceil, 
and  setting  uj®  to  be  the  average  value  of  this  function  over  the  region  inbox(t,;).  In 
Figure  3  inbox(t,_7)  would  contain  an  area- weighted  average  of  two  grid  values  while 
the  value  for  inbox(i  -f-  l,j)  is  beued  on  four  grid  values.  We  think  of  the  outbox  as 
containing  the  value  u\^'  = 

To  find  the  weights  needed  to  compute  we  must  compute  the  intersection  of  the 
inbox  with  each  neeirby  cell.  This  is  easily  accomplished  with  standard  computational 
geometry  routines.  Note  that  for  a  given  geometry  amd  grid  these  weights  need  only 
be  computed  once  at  the  beginning  of  the  computation.  They  need  not  be  recomputed 
in  each  time  step. 

Fluxes  at  other  sides.  We  now  consider  the  fluxes  F  and  G  cdong  other  sides 
of  this  cell.  These  cire  all  computed  by  similar  procedures,  so  to  be  specific  we  will 
consider  the  computation  of  the  flux  on  the  left  side  of  this  cell. 

To  compute  Fij  we  solve  two  Riemann  problems,  one  in  some  direction  (  with 
some  data  uf,  and  the  other  in  the  orthogonal  direction  7  with  data  u*,  u^.  The 
choice  of  these  directions  and  the  data  will  be  discussed  in  a  moment.  First  we  explain 
how  these  Riemann  problem  solutions  are  computed  and  used  to  define  Fij. 

Figure  4  shows  a  typical  verticed  cell  interface  and  two  orthogonal  directions  ^  and 
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Figure  4:  A  vertical  cell  interface  tind  the  and  7;-directions. 

T],  Let  9  be  the  angle  that  (  is  rotated  from  the  i-direction  (5  <  0  in  this  example). 
Suppose  we  solve  a  one-dimensional  Riemann  problem  in  the  ^-direction  with  left  and 
right  states  u^,  u*  to  obteun  the  flux  per  unit  length  per  unit  time  in  the  ^-direction. 
(To  do  this  we  rotate  the  velocity  cor  ponents  of  u|,  into  (-7]  velocity  components, 
solve  the  one- dimensioned  Riemann  problem,  and  then  rotate  the  resulting  flux  /  back 
to  x-y  velocity  components.)  Call  this  resulting  flux  f(. 

Similarly,  solving  a  one-dimensioned  Riemann  problem  in  the  7;-direction  with 
left  and  right  states  u*,  gives  the  flux  per  unit  length  per  unit  time  in  the 
7;-direction.  The  total  flux  across  the  vertical  segment  of  length  h'  is  then 

F  =  h'(f(Cos9  —  fr,sin9).  (3.19) 

This  is  the  value  we  use  for  the  flux  Fij. 

This  same  approach  has  been  used  by  others  (e.g.,  [4],  [8],  [9])  to  define  multi¬ 
dimensional  upwind  methods.  In  these  methods  the  directions  (  and  rj  are  chosen 
based  on  the  local  flow  in  an  attempt  to  use  physically  meaningful  directions  in  place 
of  the  eirtificial  coordinate  directions.  For  example,  the  direction  of  the  velocity  or 
the  pressure  gradient  might  be  used  to  define  In  our  application  we  are  only  con¬ 
sidering  cells  adjacent  to  the  boundary  and  the  relevant  directions  are  the  directions 
tangential  and  normal  to  the  wall.  We  choose  ^  to  be  the  direction  tangential  to  the 
wall  in  one  of  the  two  cells  bordering  this  interface.  Since  our  primary  concern  is  to 
maintain  stability  in  very  small  cells,  we  choose  the  smaller  of  the  two  adjacent  cells 
to  define  this  direction.  This  will  lead  to  cancellation  of  fluxes  in  tiny  cells  in  the 
same  manner  as  previously  seen  in  the  one-dimensional  example.  The  77-direction  is 
normal  to  the  ^-direction. 

Tangential  boxes.  We  must  still  specify  the  data  for  these  tangential  and  normcd 
Riemann  problems.  We  first  consider  the  tangential  problem.  Wc  use  an  approach 
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Figure  5:  (i)  Tangential  boxes  constructed  from  the  cell  interface,  (b)  Normal  boxes 
from  the  cell  interface  and  outbox(i  —  l,j). 

similar  to  the  specification  of  data  in  an  inbox  described  above.  i‘'rom  the  interface 
we  construct  boxes  that  extend  a  distance  h  in  the  ^-direction.  Figure  5a  shows  an 
example.  The  data  u*  is  obtained  by  an  area-weighted  average  of  the  values 
in  each  cell  the  box  overlaps.  In  our  current  implementation  we  assume  the  wall 
is  convex,  so  that  these  boxes  lie  entirely  within  the  computational  domeiin.  Each 
box  overlaps  at  most  two  grid  cells  and  the  weights  are  easily  calculated.  Since  the 
directions  (  and  t]  and  the  resulting  boxes  depend  only  on  the  geometry,  not  on  the 
flow  variables,  these  weights  can  again  be  calculated  once  and  for  all  as  a  preprocessing 
step. 

Normal  boxes.  Figure  5b  shows  the  normal  boxes  in  the  7;-direction.  The  box 
in  the  outward  direction  does  not  hit  the  boundary  and  overlaps  at  most  two  regular 
cells,  so  u*  is  calculated  aa  an  iirea- weighted  average  of  these  cell  values.  The  other  box 
may  extend  beyond  the  boundary.  If  so,  the  portion  lying  outside  the  computational 
domain  lies  in  one  or  more  outboxes,  the  artificial  cells  created  in  the  process  of 
computing  the  boundary  flux  Hij  described  above.  Figure  5b  shows  a  simple  example 
where  the  normal  box  intersects  only  one  cell  (i  —  l,jf)  and  outbox(x  —  l,j).  More 
generally  the  normal  box  might  intersect  two  cells  emd  their  outboxes,  as  happens 
for  example  when  we  compute  the  flux  j  which  involves  cells  (i,j)  and  (i,j  —  1). 
Moreover  the  two  outboxes  will  in  general  overlap  due  to  the  convexity  of  the  region. 
We  again  use  are2k- weighted  averaging  over  the  four  cells  in  question,  weighting  the 
values  Uij,  by  the  areas  of  intersection  and  then  dividing  by  the 

sum  of  all  these  areas. 

Cancellation.  Although  we  will  not  present  the  details  here,  it  can  be  shown 
that  this  way  of  defining  fluxes  leads  to  the  desired  cancellation  of  fluxes  in  very  smedl 
cells.  The  values  computed  at  each  of  the  three  sides  of  a  very  small  triangular 
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cell  are  nearly  the  same  because  of  our  construction.  They  differ  by  only  0{Aij)  as 
Aij  -*  0.  The  same  is  true  of  each  of  the  other  values  u*,  u"  and  so  by  Lipschitz 
continuity  of  the  fluxes  F,  G  and  H  we  obtain  the  required  cancellation.  Numerical 
results  show  stability  even  when  Aij  is  many  orders  of  magnitude  less  than  h^. 

Higher  order  methods  The  method  (3.17)  with  fltixes  (3.18)  is  only  first  order 
accurate  and  is  highly  dissipative.  In  our  previous  work  we  used  the  wave  propagation 
boundary  conditions  together  with  a  high  resolution  method  away  from  the  boundary 
and  obtained  reasonable  results  (e.g.,  [1]).  The  new  boundary  conditions  can  also 
be  applied  in  conjunction  with  a  high  resolution  method  and  gives  similar  results. 
Moreover,  with  our  new  formulation  it  appears  to  be  easier  to  improve  the  accuracy 
of  the  boundary  conditions,  allowing  us  to  obtain  higher  order  accuracy  overall.  The 
main  idea  is  to  introduce  slopes  in  each  cell  and  use  piecewise  linear  approximations 
in  place  of  piecewise  constants  to  define  the  fluxes.  Near  the  boundary  we  can  easily 
estimate  slopes  in  the  tangential  direction  along  the  wall  by  differencing  values  in  the 
inboxes  that  we  have  defined  above. 

These  improvements  are  still  being  investigated  and  will  be  reported  in  detail 
elsewhere.  Here  we  will  only  compare  results  obtained  with  the  method  as  we  have 
described  it  and  results  obtained  using  the  same  interior  method  with  the  wave  prop¬ 


agation  boundary  conditions  described  in  earlier  papers. 

4  Numerical  results 

We  show  one  representative  test  case,  a  supersonic  shock  going  around  an  expansion 
corner.  We  also  show  the  steady  state  solution  obtained  at  large  times.  The  exact 
rarefaction  wave  solution  is  a  simple  wave  and  can  be  computed  following  Section 
6.17  of  Whitham[ll],  for  example. 

The  geometry  we  use  is  the  rectangle  [0, 1.32]  x  [0,0.8]  with  a  solid  wall  at 

0.3  X  <  0.1 

y  =  0.3(1  -  (x  -  0.1)*)  0.1  <  X  <  0.7 

I  0.192  -0.36(x-0.7)  0.7  <  x  <  1.32. 

The  initial  conditions  consist  of  a  Mach  2.31  shock  at  x  =  0.06  with  left  and  right 
states 

/  =  5.1432,  =  2.04511,  uj  =  0,  =  9.04545 

and 

p"  =  1.4,  uf  =  0,  uj  =  0,  p"  =  1.0. 
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Figure  6;  Shock  propagation  results  at  t  =  0.4  (a)  Using  the  wave  propagation 
boundary  conditions,  (b)  Using  the  new  boundary  conditions. 

We  take  h  —  0.02  (66  x  40  grid)  and  a  time  step  k  =  0.002.  This  corresponds  to  a 
Courant  number  of  roughly  0.37  relative  to  the  regular  cells  with  cirea  h*.  For  the 
crude  form  of  Godunov’s  method  used  here,  the  stability  restriction  requires  Courimt 
number  less  than  0.5.  The  smallest  cells  near  the  boundary  have  an  £irea  roughly 
10-3/1*. 

Figure  6  shows  numerical  results  at  time  t  =  0.4,  bls  the  shock  is  rounding  the 
corner.  Results  obtained  with  the  wave  propagation  boimdary  conditions  are  shown 
in  Figure  6a,  while  Figure  6b  shows  the  results  obteiined  with  our  new  approach. 
These  results  are  very  similar.  Slight  discrepancies  can  be  seen  near  the  wall  just 
around  the  shock.  For  this  problem  both  sets  of  boundary  conditions  worked  well. 
We  have  also  performed  tests  on  other  problems  where  the  wave  propagation  method 
shows  instabilities  and  have  observed  no  such  difficulties  with  the  new  method. 

Figure  7  shows  the  steady  state  results  obtained  after  many  iterations  of  the  time 
dependent  code  (no  attempt  has  been  made  so  far  to  accelerate  convergence  for  steady 
state  solutions).  We  only  show  the  results  with  our  new  boundary  conditions.  The 
wave-propagation  boundary  conditions  give  very  similar  results.  We  use  a  coarser 
grid  than  in  the  previous  example  {h  =  0.4)  in  order  to  demonstrate  that  we  achieve 
reasonable  accuracy  along  the  boundary  even  with  a  relatively  coarse  piecewise  lin¬ 
ear  representation  of  the  boundary.  We  also  use  a  larger  computational  domain, 
[0,2]  X  [0, 1.6]  to  minimize  the  impact  of  the  far-field  boundaries.  The  true  solution 
is  a  rarefaction  wave  originating  from  the  portion  of  the  boundary  with  nonzero  cur¬ 
vature.  In  the  exact  solution  the  contour  lines  would  be  straight  lines.  Our  results 
are  contaminated  by  effects  from  the  far-field  boundary. 

Near  the  solid  wall  the  contour  lines  appear  to  show  a  boundary  layer.  This  is  an 
artifact  of  the  graphics  routine,  which  assumes  the  data  is  on  a  uniform  grid  at  cell 
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Figure  7:  Steady  state  results,  (a)  Pressure  contours,  (b)  Pressure  aJong  the  wall. 
The  solid  line  is  the  exact  solution,  -f  mairks  are  the  numerical  solution. 

centers.  Our  data  near  the  boundary  should  be  viewed  eis  an  approximation  to  the 
pointwise  value  at  the  center  of  mass  of  the  irregular  cell,  not  at  the  center  of  the  full 
Cartesian  cell. 

In  order  to  ex2unine  the  accuracy  at  the  wall,  Figure  7b  shows  plots  of  the  pressure 
adong  the  boundary,  plotted  against  arclength.  To  obtain  a  boundary  pressure,  the 
cell  value  Uij  and  the  reflected  value  are  used  to  solve  the  one-dimensional 

Riemann  problem  normad  to  the  wall  in  each  irregular  cell.  The  resulting  pressure  p* 
is  used  as  the  boundary  pressure.  Figure  7b  shows  these  results  and  adso  the  exact 
solution  (to  machine  precision)  calculated  using  the  theory  of  [11]. 

In  more  complicated  computations  we  use  adaptive  grid  refinement  to  obtain  high 
resolution  results  with  minimal  effort.  The  boundeury  conditions  described  here  cjm 
adso  V  5  used  in  conjunction  with  the  adaptive  Cartesian  grid  code  described  in  [1] 
and  [2]. 
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